// Mesh size
L = 400;
D = 20;
h  = 15;  // general resolution
hV = 1;   // vent resolution

// Outer box
cbl = newp; Point(cbl) = {-L/2.0,0.0,0.0,h};
cbr = newp; Point(cbr) = { L/2.0,0.0,0.0,h};
ctr = newp; Point(ctr) = { L/2.0, L ,0.0,h};
ctl = newp; Point(ctl) = {-L/2.0, L ,0.0,h};

cvl = newp; Point(cvl) = {-D/2.0,0.0,0.0,h};
cvr = newp; Point(cvr) = { D/2.0,0.0,0.0,h};

botL = newc; Line(botL) = {cbl,cvl}; 
botV = newc; Line(botV) = {cvl,cvr}; 
botR = newc; Line(botR) = {cvr,cbr}; 
rig = newc; Line(rig) = {cbr,ctr}; 
top = newc; Line(top) = {ctr,ctl}; 
lef = newc; Line(lef) = {ctl,cbl}; 

box = newreg; Line Loop(box) = {botL,botV,botR,rig,top,lef}; 

// Finally define the domain
domain = news; Plane Surface(domain) = {box};

// Vent refinement
Field[1] = MathEval;
Field[1].F = Sprintf( "min(%g,%g*(1+( (x/%g)^2 + (y/%g)^2 )))" 
                          , h,hV,       2*D,       5*D );
//Field[1].F = Sprintf( "%g*(min(abs(x)/3,1))" , h );
//Field[1].F = Sprintf( "%g/2*(0.25+1.75*(y/%g)^2)" , h , dy );
Background Field = 1;

